Load the libraries we will need for the analysis.
Load our functions.
#This function calculates PTS or PAS in the case of microsatellites.
pairwiseType <- function(x){
x <- t(x)
mm <- as(x, "dgCMatrix")
d <- tcrossprod(mm)
denom <- matrix(rep(rowSums(x), ncol(d)), ncol = ncol(d), byrow = FALSE)
denom <- denom + t(denom)
return(as.matrix(2*d/denom))
}
#This function calculates Sij and Sji (i.e., directional PTS) as in He et al 2018.
geneticSimilarity <- function(mat){
newmat <- tcrossprod(mat > 0)
newmat <- newmat/rowSums(mat > 0)
return(newmat)
}
#This function saves the PTS matrix values into a table for analysis purposes.
PTSmatrixToTable <- function(x){
diag(x) <- NA
x[upper.tri(x)] <- NA
df <- data.frame(SampleID1 = rownames(x)[row(x)], SampleID2 = colnames(x)[col(x)], PTS_score = c(x))
df <- na.omit(df)
return(df)
}
#This function saves the Genetic Similarity matrix values into a table for analysis purposes. Note that there are two GS scores for every isolate pair (i.e. directional PTS)
GSmatrixToTable <- function(x){
diag(x) <- NA
df <- data.frame(SampleID1 = rownames(x)[row(x)], SampleID2 = colnames(x)[col(x)], GS_score = c(x))
df <- na.omit(df)
return(df)
}
#This function saves the Genetic Similarity matrix values into a table for analysis purposes. Note that for the varcode analysis we only include retrospective PTS and we want to include the diagonal for the *var*code comparisons (i.e. the 1).
GSmatrixToTable_VC <- function(x){
x[upper.tri(x)] <- NA
df <- data.frame(SampleID1 = rownames(x)[row(x)], SampleID2 = colnames(x)[col(x)], GS_score = c(x))
df <- na.omit(df)
return(df)
}
ups <- fread("data/ecuador_ups_classification_final.csv", data.table = F, )
binary_all <- fread("data/ecuador_SAm_binary_matrix_final.csv", data.table = F)
binary_ecu <- fread("data/ecuador_binary_matrix_final.csv", data.table = F)
microsat_ecu <- fread("data/ecuador_binary_microsat_final.csv", data.table = F)
ecu_epi <- fread("data/ecuador_epi.csv", data.table = F)
outbreak_metadata <- read.csv("data/ecuador_varcode1_47types_metadata.csv")
matrix_all <- binary_all
rownames(matrix_all) <- matrix_all$DBLa_type
matrix_all <- matrix_all[, -1]
matrix_all <- as.matrix(matrix_all)
ecu_matrix <- binary_ecu
rownames(ecu_matrix) <- ecu_matrix$DBLa_type
ecu_matrix <- ecu_matrix[, -1]
ecu_matrix <- as.matrix(ecu_matrix)
There were a total of 195 DBLα types and 58 isolates in Ecuador.
ecu_curve <- specaccum(t(ecu_matrix), method = "rarefaction")
plot(ecu_curve, xvar = "individuals", ci.type = "polygon", xlab = "Number of DBLα sequences sampled in Ecuador", ylab = "Number of DBLα types identified in Ecuador")
#save this via Rstudio
There were a total of 543 DBLα types and 186 isolates in South America.
SAm_curve <- specaccum(t(matrix_all), method = "rarefaction")
plot(SAm_curve, xvar = "individuals", ci.type = "polygon", xlab = "Number of DBLα sequences sampled in South America", ylab = "Number of DBLα types identified in South America")
##save this via Rstudio
Because we are interested in looking at PTS or genetic similarity of parasites retrospecively, i.e. what happens after the outbreak, we will only examine pairwise comparisons of isolates retrospectively. This means that we will compare isolates “backwards” to identify how similar they are to the outbreak clone, or to anything else circulating before its identification.
genSim_ecu <- geneticSimilarity(t(ecu_matrix))
genSim_ecu[upper.tri(genSim_ecu)] <- NA
GStable_ecu <- GSmatrixToTable(genSim_ecu)
Here we will perform an analysis of the microsatellite PAS for every sample pair.
Note that for the PAS comparisons we actually don’t need the GS directional PTS for the majority of comparisons because the denominator is the same - however, we do have N=2 isolates with 6 alleles (i.e. 1 allele missing) and N=1 isolate with 5 alleles (i.e. 2 missing alleles), so will still calculate GS.
microsat_ecu <- microsat_ecu %>% filter(SampleID %in% ecu_epi$SampleID) %>%
column_to_rownames("SampleID") %>%
select(-one_of(c("MS_TA11", "MS_24901", "MS_C2M341", "MS_C3M691")))
microsat_PAS <- geneticSimilarity(microsat_ecu)
microsat_PAS[upper.tri(microsat_PAS)] <- NA
PAStable_ecuMS <- GSmatrixToTable(microsat_PAS)
We can now compare the PTS and PAS values for each pairwise comparison (retrospectively) by matching the microsat pairwise comparisons SampleID1-SampleID2 with var enabling a direct comparison.
We are interested in looking at: how correlated are predictions of recombinants by var compared to MS?
GStable_ecu_varMS <- GStable_ecu %>% left_join(PAStable_ecuMS, by = c("SampleID1", "SampleID2"))
GStable_ecu_varMS <- GStable_ecu_varMS %>% rename_at("GS_score.x", ~"PTS_score") %>%
rename_at("GS_score.y", ~"PAS_score")
correlation_plot <- GStable_ecu_varMS %>%
ggplot(aes(x = PTS_score, y = PAS_score)) +
geom_jitter(shape = 21, color = "black", size = 2) +
#stat_smooth_func(geom = "text", method = "lm", hjust = 0, parse = T) +
geom_smooth(method = "lm", se = T, color = "darkgrey") +
scale_y_continuous(breaks = seq(0, 1, by = 0.25)) +
theme_cowplot() +
ylab(expression(paste("P"["AS"], " score"))) +
xlab(expression(paste("P"["TS"], " score"))) +
background_grid(major = "xy", minor = "none")
cor.test(GStable_ecu_varMS$PTS_score, GStable_ecu_varMS$PAS_score, method = c("pearson"))
Pearson's product-moment correlation
data: GStable_ecu_varMS$PTS_score and GStable_ecu_varMS$PAS_score
t = 47.037, df = 1651, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.7353627 0.7766243
sample estimates:
cor
0.7567462
save_plot("viz/correlation_PTSvPAS.png", correlation_plot, base_width = 8, base_height = 6)
correlation_plot
There were a total of 72 (14%) out of 541 pairwise comparisons where PAS=1 but PTS ≤0.50.
There were a total of 307 (25%) out of 1232 pairwise comparisons where PAS>0.50 but PTS ≤0.50.
Now we want to examine the relationships among repertoires by visualizing the retrospective PTS values in a network.
For future reference I am building the network by following the guidelines from this tutorial) and this tutorial. I will build the network using gg syntax (i.e. similar to ggplot syntax) using the packages ggraph and tidygraph.
First we need to save our data in the proper format to build the network. The GStable_ecu becomes our edges dataset, and the epi data becomes our nodes dataset. Note: make sure to use the following code to create the edges (i.e. assigning id), if not the PTS edges between isolates may be wrong. This way we ensure that the id.x and id.y of each edge corresponds to their actual SampleID in the nodes dataset.
ecu_edges <- GStable_ecu %>% left_join(ecu_epi, by = c("SampleID1" = "SampleID")) %>% rename(id.x = id)
ecu_edges <- ecu_edges %>% left_join(ecu_epi, by = c("SampleID2" = "SampleID")) %>% rename(id.y = id)
ecu_edges <- ecu_edges %>% select(id.x, id.y, GS_score)
ecu_tidy <- tbl_graph(nodes = ecu_epi, edges = ecu_edges, directed = T)
The tibble tidy_graph object consists of both edges and nodes, either of which can be activated depending on what we need to do.
varcode_list <- ecu_epi %>% select(SampleID, varcode)
varcode_list$SampleID <- as.factor(varcode_list$SampleID)
varcode_list$varcode <- as.factor(varcode_list$varcode)
varcode_list <- varcode_list %>% column_to_rownames("SampleID")
varcode_colors <- list(varcode = c("varcode1" = "#fb9a99", "varcode2" = "#e538a9","varcode3" = "#48B24F", "varcode4" = "#30d5c8", "varcode5" = "#E4B031", "varcode6" = "#3090d5", "varcode7" = "#CAD93F", "varcode8" = "#9ac4b3", "varcode9" = "#a880bb"))
ecu_varcodes_network <- ecu_tidy %>% activate(edges) %>% filter(GS_score >= 0.9) %>%
ggraph(layout = "fr") +
geom_edge_link(color = "darkgrey", alpha = 0.5) +
geom_node_point(aes(fill = varcode), shape = 21, color = "black", size = 4) +
#geom_node_text(aes(label = SampleID), repel = TRUE) +
labs(edge_width = "SampleID",
fill = "") +
scale_fill_manual(values = varcode_colors$varcode, labels = c("varcode1 (N=36)", "varcode2 (N=3)", "varcode3 (N=3)", "varcode4 (N=1)", "varcode5 (N=1)", "varcode6 (N=3)", "varcode7 (N=8)", "varcode8 (N=2)", "varcode9 (N=1)")) +
theme_graph()
save_plot("viz/ecuador_varcode_network.png", ecu_varcodes_network, base_width = 8, base_height = 6)
ecu_varcodes_network
| varcode | n |
|---|---|
| varcode1 | 36 |
| varcode2 | 3 |
| varcode3 | 3 |
| varcode4 | 1 |
| varcode5 | 1 |
| varcode6 | 3 |
| varcode7 | 8 |
| varcode8 | 2 |
| varcode9 | 1 |
timepoints <- c("2013", "2014", "2015")
totalnum_varcodes <- c("3", "8", "4")
varcodes_bytime <- data.frame(timepoints, totalnum_varcodes)
The number of varcodes identified in each timepoint was as follows:
| timepoints | totalnum_varcodes |
|---|---|
| 2013 | 3 |
| 2014 | 8 |
| 2015 | 4 |
.
varcode <- c("varcode1", "varcode2", "varcode3", "varcode6", "varcode7", "varcode8", "varcode4", "varcode5", "varcode9")
vc2013 <- c(30, 2, 1, NA, NA, NA, NA, NA, NA)
vc2014 <- c(4, 1, 2, 2, 6, 2, 1, 1, NA)
vc2015 <- c(2, NA, NA, 1, 2, NA, NA, NA, 1)
varcodes_temporal <- data.frame(varcode, vc2013, vc2014, vc2015)
varcodes_temporal <- varcodes_temporal %>% rename_at("vc2013", ~"2013") %>% rename_at("vc2014", ~"2014") %>% rename_at("vc2015", ~"2015")
varcodes_timeline <- varcodes_temporal %>%
melt(id.var = "varcode") %>%
rename_at("variable", ~"Year") %>%
rename_at("value", ~"Frequency") %>%
mutate(varcode = factor(varcode, levels = c("varcode9", "varcode8", "varcode7", "varcode6", "varcode5", "varcode4", "varcode3", "varcode2", "varcode1"))) %>%
ggplot(aes(Year, factor(varcode))) +
geom_point(aes(size = Frequency, color = varcode)) +
geom_text(aes(label = Frequency), size = 2.5, vjust = 1.5, hjust = -1) +
scale_size_continuous(name = "Number of Pf isolates",
breaks = c(1, 10, 30),
labels = c("1", "10", "30")) +
scale_color_manual(values = varcode_colors$varcode, labels = varcode_list) +
annotate("segment", x = 1, xend = 3, y = 9, yend = 9, colour = "#fb9a99") +
annotate("segment", x = 1, xend = 2, y = 8, yend = 8, colour = "#e538a9") +
annotate("segment", x = 1, xend = 2, y = 7, yend = 7, colour = "#48B24F") +
annotate("segment", x = 2, xend = 3, y = 4, yend = 4, colour = "#3090d5") +
annotate("segment", x = 2, xend = 3, y = 3, yend = 3, colour = "#CAD93F") +
ylab("") +
theme_bw()
save_plot("viz/varcodes_timeline.png", varcodes_timeline)
varcodes_timeline
Now we can calculate the average days a varcode persisted.
ecu_epi$DateCollected <- as.Date(ecu_epi$DateCollected)
#ecu_epi %>% filter(SampleID == "ECPf003" | SampleID == "ECPf076") %>% select(varcode, DateCollected)
vc1_duration <- data.frame(varcode = "varcode1", FirstID = subset(ecu_epi, SampleID == "ECPf003")$DateCollected, LastID = subset(ecu_epi, SampleID == "ECPf076")$DateCollected)
vc2_duration <- data.frame(varcode = "varcode2", FirstID = subset(ecu_epi, SampleID == "ECPf024")$DateCollected, LastID = subset(ecu_epi, SampleID == "ECPf045")$DateCollected)
vc3_duration <- data.frame(varcode = "varcode3", FirstID = subset(ecu_epi, SampleID == "ECPf035")$DateCollected, LastID = subset(ecu_epi, SampleID == "ECPf053")$DateCollected)
vc6_duration <- data.frame(varcode = "varcode6", FirstID = subset(ecu_epi, SampleID == "ECPf051")$DateCollected, LastID = subset(ecu_epi, SampleID == "ECPf072")$DateCollected)
vc7_duration <- data.frame(varcode = "varcode7", FirstID = subset(ecu_epi, SampleID == "ECPf056")$DateCollected, LastID = subset(ecu_epi, SampleID == "ECPf067")$DateCollected)
vc_duration <- bind_rows(vc1_duration, vc2_duration, vc3_duration, vc6_duration, vc7_duration)
vc_duration <- vc_duration %>% mutate(duration_days = difftime(LastID, FirstID, units = c("days"))) %>%
mutate(duration_yr = round(duration_days/365, 2)) %>%
mutate(duration_months = round(duration_yr*12, 2))
The varcodes persisted for the following amount of time:
| varcode | FirstID | LastID | duration_days | duration_yr | duration_months |
|---|---|---|---|---|---|
| varcode1 | 2013-01-28 | 2015-05-01 | 823 days | 2.25 days | 27.00 days |
| varcode2 | 2013-07-15 | 2014-01-21 | 190 days | 0.52 days | 6.24 days |
| varcode3 | 2013-10-24 | 2014-05-28 | 216 days | 0.59 days | 7.08 days |
| varcode6 | 2014-03-13 | 2015-11-08 | 605 days | 1.66 days | 19.92 days |
| varcode7 | 2014-07-08 | 2015-01-19 | 195 days | 0.53 days | 6.36 days |
| min_days | med_days | avg_days | max_days |
|---|---|---|---|
| 190 days | 216 days | 405.8 days | 823 days |
| min_months | med_months | avg_months | max_months |
|---|---|---|---|
| 6.24 days | 7.08 days | 13.32 days | 27 days |
| min_yr | med_yr | avg_yr | max_yr |
|---|---|---|---|
| 0.52 days | 0.59 days | 1.11 days | 2.25 days |
ecu_tidy %>% activate(nodes) %>% filter(LocationCode == "San Lorenzo, Esmeraldas") %>%
activate(edges) %>% filter(GS_score >= 0.9) %>%
ggraph(layout = "fr") +
geom_edge_link(color = "darkgrey", alpha = 0.5) +
geom_node_point(aes(fill = varcode), shape = 21, color = "black", size = 4) +
#geom_node_text(aes(label = SampleID), repel = TRUE) +
labs(edge_width = "SampleID",
fill = "") +
scale_fill_manual(values = varcode_colors$varcode) +
theme_graph()
ecuador_recombinants_network <- ecu_tidy %>% activate(edges) %>% filter(GS_score >= 0.5) %>%
ggraph(layout = "fr") +
geom_edge_link(color = "darkgrey", alpha = 0.5) +
geom_node_point(aes(fill = varcode), shape = 21, color = "black", size = 4) +
scale_fill_manual(values = varcode_colors$varcode, labels = varcode_list) +
#geom_node_text(aes(label = SampleID), repel = TRUE) +
labs(edge_width = "SampleID") +
theme_graph()
save_plot("viz/ecuador_recombinants_network.png", ecuador_recombinants_network, base_width = 8, base_height = 6)
ecuador_recombinants_network
ecuMSedges <- PAStable_ecuMS %>% left_join(ecu_epi, by = c("SampleID1" = "SampleID")) %>% rename(id.x = id)
ecuMSedges <- ecuMSedges %>% left_join(ecu_epi, by = c("SampleID2" = "SampleID")) %>% rename(id.y = id)
ecuMSedges <- ecuMSedges %>% select(id.x, id.y, GS_score)
ecuMS_tidy <- tbl_graph(nodes = ecu_epi, edges = ecuMSedges, directed = T)
ecuador_MSclusters_network <- ecuMS_tidy %>% activate(edges) %>% filter(GS_score >= 0.90) %>%
ggraph(layout = "fr") +
geom_edge_link(color = "darkgrey", alpha = 0.5) +
geom_node_point(aes(fill = varcode), shape = 21, color = "black", size = 4) +
scale_fill_manual(values = varcode_colors$varcode, labels = varcode_list) +
#geom_node_text(aes(label = SampleID), repel = TRUE) +
labs(edge_width = "SampleID") +
theme_graph()
save_plot("viz/ecuador_MSclusters_network.png", ecuador_MSclusters_network, base_width = 8, base_height = 6)
ecuador_MSclusters_network
Allowing for 1 allele diff
ecuador_MSclusters_network80 <- ecuMS_tidy %>% activate(edges) %>% filter(GS_score >= 0.80) %>%
ggraph(layout = "fr") +
geom_edge_link(color = "darkgrey", alpha = 0.5) +
geom_node_point(aes(fill = varcode), shape = 21, color = "black", size = 4) +
scale_fill_manual(values = varcode_colors$varcode, labels = varcode_list) +
#geom_node_text(aes(label = SampleID), repel = TRUE) +
labs(edge_width = "SampleID") +
theme_graph()
save_plot("viz/ecuador_MSclusters_network_PTS80.png", ecuador_MSclusters_network80, base_width = 8, base_height = 6)
ecuador_MSclusters_network80
ecuador_MSclusters_supp <- ecuador_MSclusters_network + ecuador_MSclusters_network80 +
plot_layout(guides = "collect") + plot_annotation(tag_levels = "a")
save_plot("viz/ecuador_MSclusters_networks_both.png", ecuador_MSclusters_supp, base_width = 12, base_height = 6)
ecuador_MSclusters_supp
ecuador_MSrecomb_network <- ecuMS_tidy %>% activate(edges) %>% filter(GS_score >= 0.5) %>%
ggraph(layout = "fr") +
geom_edge_link(color = "darkgrey", alpha = 0.5) +
geom_node_point(aes(fill = varcode), shape = 21, color = "black", size = 4.5) +
scale_fill_manual(values = varcode_colors$varcode, labels = varcode_list) +
#geom_node_text(aes(label = SampleID), repel = TRUE) +
labs(edge_width = "SampleID") +
theme_graph()
save_plot("viz/ecuador_MSrecomb_network.png", ecuador_MSrecomb_network, base_width = 8, base_height = 6)
ecuador_MSrecomb_network
This very clearly shows that describing transmission dynamics using microsatellites is not sufficient to identify recent genetic exchanges and recent recombination/transmission events since everything is connected by MS. Since var genes are rapidly evovling compared to MS, they allow us to “track” recombination and transmission dynamics in space and time.
###Trying to order the clusters by "year"
pheatmap(t(ecu_matrix), col = c("white", "black"),
annotation_row = varcode_list,
annotation_colors = varcode_colors,
fontsize_row = 5, fontsize_col = 5,
cluster_rows = T, cluster_cols = F,
legend = F,
treeheight_col = 0,
show_colnames = F,
show_rownames = F,
#filename = "viz/varcode_heatmap.png",
width = 8,
height = 4)
###Trying to order the clusters by "year"
pheatmap(microsat_ecu, col = c("white", "black"),
annotation_row = varcode_list,
annotation_colors = varcode_colors,
fontsize_row = 5, fontsize_col = 5,
cluster_rows = T, cluster_cols = F,
legend = F,
treeheight_col = 0,
show_colnames = F,
show_rownames = F,
#filename = "viz/microsat_heatmap.png",
width = 8,
height = 4)
register_google(key = "INSERT_KEY")
geocode("Ecuador")
map <- get_googlemap(center = "Ecuador", zoom = 7, maptype = "terrain", style = "feature:poi.business|visibility:off&style=feature:road|element:labels.icon|visibility:off&style=feature:transit|visibility:off")
ecu_map <- ggmap(map) + xlab("Longitude") + ylab("Latitude")
ecu_map <- ecu_map + scalebar(x.min = -77, x.max = -75,
y.min = -4.8, y.max = -4.7,
dist = 100,
dist_unit = "km",
transform = T,
model = "WGS84",
#st.bottom = F,
st.dist = 1,
height = 1,
st.size = 3) +
annotation_north_arrow(location = "br",
pad_x = unit(0.4, "in"),
pad_y = unit(0.4, "in"),
style = north_arrow_fancy_orienteering)
ecu_map_zoom <- ggmap(map) + xlab("Longitude") +
ylab("Latitude") +
ylim(-2.25, 1.4) +
scalebar(x.min = -77, x.max = -75,
y.min = -2.1, y.max = -2.0,
dist = 100,
dist_unit = "km",
transform = T,
model = "WGS84",
#st.bottom = F,
st.dist = 1,
height = 1,
st.size = 3) +
annotation_north_arrow(location = "br",
pad_x = unit(0.9, "in"),
pad_y = unit(0.4, "in"),
style = north_arrow_fancy_orienteering)
ecu_map_piecharts <- ggmap(map) + xlab("Longitude") +
ylab("Latitude") +
ylim(-2.25, 1.6) +
scalebar(x.min = -77, x.max = -75,
y.min = -2.1, y.max = -2.0,
dist = 100,
dist_unit = "km",
transform = T,
model = "WGS84",
#st.bottom = F,
st.dist = 1,
height = 1,
st.size = 3) +
annotation_north_arrow(location = "br",
pad_x = unit(0.1, "in"),
pad_y = unit(3.2, "in"),
style = north_arrow_fancy_orienteering)
ecuador_gmap_locations <- ecu_map + geom_point(data = ecu_epi, aes(x = lon, y = lat, fill = LocationCode), shape = 21, color = "black", size = 3) +
scale_fill_manual(values = loc_cols,
labels = locationLabels$LocationCode) +
labs(fill = "Location")
save_plot("viz/ecuador_googlemap_samplinglocations.png", ecuador_gmap_locations)
ecuador_gmap_locations
ecu_epi <- ecu_epi %>%
mutate(Year = case_when(DateCollected < "2013-12-31" ~"2013",
DateCollected < "2014-12-31" ~"2014",
DateCollected < "2015-12-31" ~"2015",
TRUE ~"CHECK"))
barplot_samples_perlocation <- ecu_epi %>%
group_by(LocationCode, Year) %>%
tally() %>%
ggplot(aes(x = Year, y = n, fill = LocationCode)) +
geom_bar(stat = "identity") +
scale_fill_manual(values = loc_cols,
labels = locationLabels$LocationCode) +
background_grid(major = "xy", minor = "none") +
labs(x = "Year",
y = "Number of samples",
fill = "Location") +
theme(legend.position = "none") +
theme_cowplot() +
background_grid(major = "xy")
save_plot("viz/ecuador_barplot_sampling.png", barplot_samples_perlocation)
barplot_samples_perlocation
seg_data <- ecu_edges %>% filter(GS_score >= 0.5)
seg_data <- seg_data %>% left_join(ecu_epi, by = c("id.x"="id"))
seg_data <- seg_data %>% left_join(ecu_epi, by = c("id.y"="id"))
seg_data <- seg_data %>% group_by(varcode.x, varcode.y, lon.x, lat.x, lon.y, lat.y) %>% tally()
ecu_map_zoom +
geom_segment(data = seg_data,
aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = 0.1, size = n),
linejoin = "mitre") +
geom_point(data = ecu_epi,
aes(x = lon, y = lat, fill = varcode),
shape = 21, color = "black", size = 3.5,
position = position_jitter(h = 0.1, w = 0.1)) +
scale_size_area() +
scale_fill_manual(values = varcode_colors$varcode,
labels = varcode_list)
seg_data1 <- ecu_edges %>% filter(GS_score >= 0.5)
seg_data1 <- seg_data1 %>% left_join(ecu_epi, by = c("id.x"="id"))
seg_data1 <- seg_data1 %>% left_join(ecu_epi, by = c("id.y"="id"))
seg_data1 <- seg_data1 %>% rename_at("DateCollected.x", ~"DateCollected")
googlemapSamplesRTNet <- ecu_map_zoom +
geom_segment(data = seg_data1,
aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = 0.1, size = 0.5),
linejoin = "bevel") +
geom_point(data = ecu_epi,
aes(x = lon, y = lat, fill = varcode),
shape = 21, color = "black", size = 3.5,
position = position_jitter(h = 0.1, w = 0.1)) +
scale_size_area() +
scale_fill_manual(values = varcode_colors$varcode,
labels = varcode_list) +
scale_alpha(guide = "none") +
scale_size(guide = "none") +
transition_states(DateCollected,
transition_length = 1,
state_length = 5) +
labs(title = "",
subtitle = 'Samples collected on {closest_state}') +
shadow_mark(size = 3)
googlemapSamplesRTNet_anim <- animate(googlemapSamplesRTNet, width = 8, height = 8, res = 500, units = "in", duration = 15, fps = 2, detail = 5, render = gifski_renderer())
anim_save("viz/spatial_network_realtime_googlemap.gif", googlemapSamplesRTNet_anim)
googlemapSamplesRTNet_anim <- animate(googlemapSamplesRTNet, width = 8, height = 8, res = 250, units = "in", duration = 15, fps = 2, detail = 5, render = gifski_renderer())
anim_save("viz/spatial_network_realtime_googlemap_lowres.gif", googlemapSamplesRTNet_anim)